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Abstract 

Background: Glomerulus detection is a key step in histopathological evaluation of microscopy images of 

kidneys. However, the task of automatic detection of glomeruli poses challenges due to the disparity in sizes and 
shapes of glomeruli in renal sections. Moreover, extensive variations of their intensities due to heterogeneity in 
immunohistochemistry staining are also encountered. 

Despite being widely recognized as a powerful descriptor for general object detection, the rectangular histogram of 
oriented gradients (Rectangular HOG) suffers from many false positives due to the aforementioned difficulties in 
the context of glomerulus detection. 

Results: A new descriptor referred to as Segmental HOG is developed to perform a comprehensive detection 
of hundreds of glomeruli in images of whole kidney sections. The new descriptor possesses ffexible blocks that 
can be adaptively fitted to input images to acquire robustness to deformations of glomeruli. Moreover, the novel 
segmentation technique employed herewith generates high quality segmentation outputs and the algorithm is assured 
to converge to an optimal solution. Gonsequently, experiments using real world image data reveal that Segmental 
HOG achieves significant improvements in detection performance compared to Rectangular HOG. 

Conclusion: The proposed descriptor and method for glomeruli detection present promising results and is 

expected to be useful in pathological evaluation. 


1 Background 

Renal glomeruli provide a filtration barrier that retains higher molecular weight proteins in blood circulation. In 
various renal diseases, damage of the glomerular filtration barrier can be observed as protein leakage into urine, 
known as proteinuria. Therefore, the pathological changes in renal glomeruli of animal disease models can provide 
important information in screening compounds that target such diseases. 

Our goal is to perform high-throughput detection of glomeruli in huge microscopy images of animal disease models, 
whose sizes run up to the order of 10^ pixels. While existing studies about automatic analysis of glomeruli in microscopy 
images of kidneys are present ESI US], target images in these works are from human biopsy samples with relatively 
small sizes, and are not suitable for our purpose. 

Compared to general object detection tasks, there are two particular obstacles in the case of glomerulus detection. 
The first obstacle arises from the non-rigid sizes and shapes of the targets within the images. Indeed, the sizes of 
glomeruli are stable in vivo^ although they swell in unfavorable situations, e.g. hypertension [6] and diabetes EH , to 
some degree. Also, the sizes of glomeruli in a whole-kidney-section image could vary depending on which part of a 
glomerulus the cross section passes through. The shapes of the glomeruli are almost spherical, making the boundaries 
circular. To obtain the boundaries, one might try to fit an ellipse to each glomerulus. However, this approach yields 
large estimation errors because each glomerulus is deformed to some extent. 

The second difficulty arising in glomerulus detection task is the high variation of the intensities. In histological 
evaluation, immunohistochemistry is usually used to demonstrate the distribution and location of proteins in sections. 
In our target images, sections are immunostained for desmin, a known glomerular injury marker. As a result, some 
glomeruli are stained, and some are not. Since many glomeruli are partly stained, yielding heterogeneously stained 
glomeruli, detection is more complicated. Furthermore, the stained tissues in the kidneys are not only from glomeruli 
but also from other tissues such as blood vessels. 

To check the existence of a glomerulus at each location in a whole-kidney-section image, the sliding window 
technique [151 m [ISl EO] is employed. Using this procedure, a frame goes over the input image to check at every 
possible location whether the target object exists, then, a descriptor of the sub-image is extracted. 
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Rectangular HOG (R-HOG) [4], a widely used and recognized efficient descriptor for object detection in the field 
of computer vision, is a potentially suitable candidate descriptor for glomeruli. It has the capacity to capture the 
information of the magnitudes of the gradients in the image. Therefore, it is robust to the change in intensities caused 
by the heterogeneity of the stained levels. Glomeruli are known to be composed of tightly packed cells, resulting to 
high gradients in images. Thus, a natural approach would be to use the magnitudes of these gradients as features for 
the glomeruli. While we have also previously attempted to directly exploit this attribute, we have found the detection 
performance to be poor, resulting to many false positives and low recall. In addition to the magnitudes of gradients, 
their directions are also important to distinguish glomeruli from the other tissues. Using R-HOG descriptors obtained 
from both the magnitudes and the directions of the gradients, glomeruli detection performance results to recall values 
high enough to be useful for pathological evaluation. However, it appeared that R-HOG still suffers from a considerable 
amount of false positives [30151. 

The high number of false positives from the previous studies I30i can be ascribed to the condition that the 
standard HOG such as the R-HOG has a rigid block division. Due to this rigidity, there are instances when a block is 
inside the glomerular area, and outside on another. Thus, extracted features from each block contain large variances, 
and robustness to the deformations of glomeruli is lost. Although there are several other known local descriptors 
such as scale-invariant feature transform (SIFT) image features [12], Haar-like features [20|, and local binary patterns 
(LBP) H], these do not possess a solution to be robust to the deformation of glomeruli for similar reasons. 

In this study, we introduce a flexible block division to the HOG descriptor to improve the detection performance and 
reduce the number of false positives. A new feature, which we refer to as the Segmental HOG (S-HOG) descriptor, is 
proposed for glomerulus detection. The block division of S-HOG is based on the estimated boundary of the glomerulus 
that is obtained via a segmentation algorithm, which is also developed in this work. This renders the division of blocks 
to be more adaptable than the rigid block division of R-HOG, and allows feature vectors to clearly differentiate between 
the inside and the outside of the glomerulus. Moreover, since blocks are always within the glomerular area, gradient 
information in the same block between two glomeruli is expected to be more similar. Experiments conducted reveal 
that the number of false positives was halved, keeping almost all true positives when using S-HOG compared to the 
R-HOG. 


Related Works 

Segmentation is an important step to extract the S-HOG descriptors. Recent works on segmentation of glomeruli has 
been sparse mini. Nevertheless, there has been some research regarding segmentation of specific organs in general 
biomedical images, including region growing [7], level set method [3|, and active contour model [Tnil2Ill23- Majority of 
them are semi-automatic and require users’ intervention, possess no guarantee of optimality [EliniES, and are highly 
dependent on the initial solution provided by users as input. On the other hand, the segmentation algorithm developed 
in this study is ensured theoretically to obtain the optimal solution, producing high quality segmentation. In addition 
to the above-mentioned, most recent attempts include using deep learning [2|. Deep learning typically requires great 
computational and time resources, whereas the proposed algorithm can work even on a standard personal computer 
or a laptop. 

The algorithm developed by Kvarnstrom et al. m is relevant to the proposed segmentation technique. Their 
algorithm for cell contour recognition is based on a dynamic program, where they first estimated the cell centers and 
constructed a ray from the center to each of m directions, where m = 32. Then they computed the boundary likeliness 
at n points on each ray, where they set n = 30. Their algorithm finds a smooth contour by taking a point on each ray 
to connect them. To get a closed contour, they presented two algorithms. The first algorithm poses n sub-problems 
where in each sub-problem, the initial point and the endpoint are the same. We shall refer to this algorithm here 
as the exhaustive dynamie program (EDP). Their second algorithm is a heuristic method that is faster than the first 
scheme, but possesses no guarantee for global optimality. In this study, we developed a new segmentation algorithm 
which we will refer to as divide & eonquer dynamie program (DGDP). Gompared to Kvarnstrom et al’s algorithms [TT] . 
the advantages of the DGDP algorithm are as follows: 

• DGDP is much faster than EDP, yet an exact optimal solution is always obtained. 

• The boundary likeliness function is trained with a machine learning technique to precisely estimate boundaries 
of glomeruli. 

2 Methods 

In the proposed method, a new descriptor, S-HOG, is introduced to detect glomeruli in kidney microscopic images. 
Segmentation of glomeruli is needed to extract the S-HOG descriptor. Eor fast exhaustive detection of glomeruli from 
large microscopic images, pre-screening is performed with R-HOG which does not require prior segmentation. The 
proposed method consists of the three stages (EigureHj): 
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• pre-screening stage, 

• segmentation stage, and 

• classification stage. 

In each stage, a support vector machine (SVM) [181 US] is used with a different type of HOG descriptor, resulting in 
three SVMs in total. To obtain the S-HOG descriptor, we perform segmentation of glomeruli from the sub-images 
that passed the pre-screening (Figured]). Hereinafter, we present the details of each stage, and discuss how training 
datasets for each SVM are constructed and the materials used in the experiments at the end of this section. 

Pre-Screening 

In the pre-screening stage, candidate glomeruli are detected from a kidney microscopy image using the sliding window 
technique. The window size is set to 200 x 200 in our experiments. R-HOG features, which are 512-dimensional 
vectors based on our selected parameter values, are extracted and judged by SVM, and non-maximal suppression is 
then performed to obtain candidate glomeruli. 

Segmentation 

Segmentation of glomeruli is performed on sub-images that passed the pre-screening. In the segmentation algorithm, 
the boundary of a glomerulus is represented by an m-sided polygon whose m vertices are restricted to he on m line 
segments, respectively. The m line segments are placed uniformly, as outlined by the dotted lines in Figure [2)3 for the 
case where m = 36. To determine the location of the vertex on each line segment, the sliding window technique is 
employed again. 0 The window sweeps through the line segment and computes the boundary likeliness at n locations 
on the line segments. In Figure [2)3, the boundary likeliness Li {i = 1,... ,m) is computed at every dotted location. 
How Li is computed is discussed at the end of this Segmentation subsection. We set the length of the line segment to 
63 pixels, where the endpoint closest to the center of the image is 17 pixels away from the center. The length between 
adjacent dots along a line segment is equal to 3 pixels, resulting to n = 22 dots on each line segment. In total, the 
boundary likeliness is computed at mn{= 36 x 22 = 792) locations. To determine the vertices of the m-sided polygon, 
one might consider naively locating the points that achieve the highest boundary likeliness on each line segment. 
However, this approach often yields an extremely zigzag boundary. 

To obtain a smoother boundary, we impose a constraint that suppresses distant adjacent vertices. We then establish 
the following maximization problem: 

m 

max 

subj to \pi-p 2 \<^, 

where L* : {1 ,... ,n} —> M denotes the boundary likeliness function obtained by the sliding window technique, and 
Pi e , n} {i = 1 ,..., m) is a location on the i-th line segment, where the line segment is discretized into n points 

numbered with a natural number. For instance, when = 1, the Tth vertex is at the endpoint of the i-th line segment 
closest to the center, and the vertex can move from this endpoint to the other endpoint with increasing value of Pi. 
Li{p) is the boundary likeliness at p-th location in the i-th line segment. The location pi on the i-th line segment 
is more likely to be the boundary with larger Li{pi) value. In our experiments, we set ^ = 1. We shall denote the 
objective function by J(p), i.e., J{p) := Li{pi). 

Kvarnstrom et al. tackled a similar optimization problem, and proposed an algorithm which we called EDP. In this 
study, a faster algorithm named DCDP is developed. The two algorithms, EDP and DCDP, are detailed as follows. 

Exhaustive Dynamic Program (EDP). Due to the last constraint \pm —pi\ < ^5 the standard dynamic program 
cannot directly solve optimization problem ([T]). To make the problem tractable, EDP divides problem ([T]) into n 
sub-problems, where the value of Pm is fixed to one of {1 ,... ,n} in each sub-problem. Once all the n sub-problems 
are solved, the optimal solution of the original problem m is obtained by taking the solution that yields the largest 
objective value among n sub-problems. Each sub-problem can be solved by a dynamic program that takes 0{nm(;) 
computation. This algorithm computes an optimal solution of the k-th. sub-problem, which is equivalent to the 
original problem ([T]) to which the constraint Pm = k is added. Hence, the maximization problem ([T]) can be solved 
with computation. 

The new algorithm DCDP also takes 0{n^m<;) computational time in worst case, although the new algorithm 
solves the same problem much faster than EDP, as presented in the Results section. 

^Sliding windows are used in both pre-screening and segmentation. 


wrt pi,... e { 1 ,... ,n}, 
■■,\Pm-l -Pm\ < \Pm “ Pi I < 
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Divide Conquer Dynamic Program (DCDP). Observe that problem m can be solved in 0{nm(^) compu¬ 
tational time by a dynamic program if one of the constraints \pm — Pi\ < ^ is disregarded. The idea to devise the new 
algorithm is based on the following fact: Suppose is an optimal solution that maximizes J(-) without the constraint 
\Pm — Pi\ < Then if pg is a feasible solution for the original problem ([1]), it is also an optimal solution. 

To express this idea mathematically, let us define 

S{I) :={p = {pi,... ,pm) e N™ I 

Pm e X, \pi -P2\ <?,•••, \Pm-l -Pm\ < \Pm “ Pi | < ?}• 

for X C Nn, where Nn := {1 ,... ,n}. Note that 5(Nn) is the feasible region of the original problem ([T]). The goal of 
DCDP is to find an optimal solution 


p^ G argmax J (p ). 

peSiNr^) 

Dynamic program (DP) cannot solve this problem directly due to the existence of the constraint \pm —pi\<^- To 
use DP, we consider finding the maximizer of J{p) from a relaxed region Si,{Nn) where 5 l(X) is defined as 

5l(X) :={p={p,,...,pm)GN^\ 

Pm e X, Pi e X + \pi -P2I < ?,•••, |Pm-l -Pm\ < ?}, 

where the operator + denotes that for any two sets I and J^ X^ J \= {i ^ j \ i ^ j G J}. Note that S{X) C 5l(X). 
The strategy of DCDP is to first find the solution of the relaxed problem, 

Po ^ argmax J(p) 
pe5L(N^) 

and then check the feasibility: if pg G 5(Nn), then pg is the optimal solution of the original problem ([T]). If Pg ^ tS(Nn), 
then the set is divided into Xi and X 2 (i.e. Xi UX2 = N^), and the following two sub-problems are solved: 

p\ G argmax J(p) and p\ G argmax J(p). 
pes{Xi) pes{X2) 

Notice that the original feasible region tS(Nn) is the sum of the two regions, S{Xi) and S{X 2 ). Therefore, we can take 
either of the two solutions, p\ and P2, which has the larger objective value. DCDP employs a divide and conquer 
approach that repeatedly applies the above strategy to sub-problems. The basic approach of DCDP is summarized in 
Algorithm [H Invoking the function DCDP _B ASIC (Nn) yields the optimal solution of the original problem. Here, the 
function (Xi,X2 ) := Split(Xo) divides the set Xq into two exclusive non-empty subsets, Xi and X 2 . 


Algorithm 1 po = DCDP .Basic (Xq). 

Input: Xo C Nn- 

Output: Po e argmaxpg 5 ( 2 :^) J(p). 

1: Po e argmaxpe5^(2:^) J(p); 

2: if Pg G tS(Xg), then return; 

3: (Xi,X2) := SPLlT(Xg); 

4: Pi := DCDP_Basic(Xi); 

5: P2 := DCDP_Basic(X2); 

6: G argmax^^ii 2} 

7: Po :=PiN 


The first step po G argmax^^^^^^j^^ J(p) can be performed in 0{nm<;) computational time. An instance of the 
dynamic program is given in Algorithmic Note that po G tS(Xg) is always ensured if the cardinality of Xg is one since 
the relaxed region is reduced to the unrelaxed region (i.e. S{{h}) = Si,{{h})). The function DCDP .Basic is invoked, 
at most, (2n — 1) times. This implies that the computational time in worst case is As will be shown in 

the Results section, we empirically found that the number of invoking the function recursively is much smaller than 
(2n- 1). 

To accelerate the DCDP algorithm, pruning steps are added. For the implementation of SPLlT(Xg), we considered 
three schemes: Half Split, Max Split, and Adap Split. The pruning steps and the three schemes of Split (Xg) are 
detailed in Sections El and m as well as the resulting accelerated DCDP algorithm. A mathematical proof using case 
analysis showing that our algorithm always obtains an optimal solution is also given in Section O 
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Algorithm 2 0{nm(;) Dynamic Program for J{p) 

Input: X C Nn- 

Output: p G argmaxpg 5 ^(j) J(p) 

1: Initialize all entries in the n x m matrix Q with — oo. 

2: for j el + +?] n N„ do 

3: Q{j,l) := Li{j); 

4: end for 

5: for t = 2, ..., m do 

6: if t < m, then Xf := N^else Xf : = X; 

7: for i G Xt do 

8: j* := argmax{(3(j, t - 1) \ j e [i - <;,i + <i] n N„}; 

9: Q(i,t) := Lt(i)+Q(i*,t-1); P{i,t) := j^,; 

10: end for 

11: end for 

12: e argmaxjgj( 3 (i,m); 

13: t := m — 1; 

14: while t > 1 do 

15: := P(p^+i,t + 1); 

16: end while 


Computation of Boundary Likeliness The sliding window technique is employed in order to determine 

the vertices of the m-sided polygon described above. The window size in this stage is set to 30 x 15 pixels, and the 
windows sweep through the line segments (Figure [2)3). Each time the sliding window moves, a feature descriptor is 
computed from the window and is applied with Linear SVM to compute the SVM score, which is what we refer to as 
the boundary likeliness Li{'). The SVM scores of m vertices, 1/^(1),..., Li{n), are then obtained for i = 1,..., m, and 
integrated in the maximization problem ©■ 

The HOG feature is adopted as the descriptor to compute the boundary likeliness. Each window is divided into 
three blocks as shown in Eigure[2)3. This division design is from an observation that some glomeruli are surrounded 
with a thick Bowman’s capsule, and the middle block is expected to capture this glomerular capsule. The statistics of 
nine discretely oriented gradients are computed in each block, producing a 27-dimensional feature vector. 

Classification with S-HOG descriptor 

Candidate glomeruli obtained via pre-screening are classified using the proposed S-HOG descriptor. S-HOG exploits 
the glomerulus boundary located in the segmentation stage to generate 24 non-overlapping blocks, as shown in Eig- 
ure[3fc. 

Various types of glomeruli are contained in kidney microscopy images, some of them surrounded with a thick 
Bowman’s capsule. To effectively exploit this characteristic, the circle containing a candidate glomerulus is divided 
into three zones: Inner zone, middle zone, and outer zone. We divide the circle into eight disjoint sectors, and take 
the intersection of each zone and each sector to get the 24 non-overlapping blocks (Eigure[3b), and gradients are then 
histogrammed for each block (Eigure [3li). In our experiments, we employed nine discretized oriented gradients, and 
SVM is applied to S-HOG feature vectors to discriminate between glomeruli and non-glomeruli. 

Construction of Training Data 

A total of three linear SVMs are used, one for the pre-screening, segmentation, and classification stage, respectively. 
A training dataset is required for each of the three SVMs. Details on the construction of each training data set are 
given subsequently. 

Training Data for Pre-Screening Stage. Each example in the training data for pre-screening is a 200 x 200 
sub-image. A positive example contains a glomerulus in the center of the sub-image, while a negative example does 
not. To gather these samples, the locations of glomeruli within the whole-kidney-section images used for training are 
first annotated manually. Small glomeruli whose diameters are under 50 pixels were ignored. Positive examples are 
the sub-images from 200 x 200 bounding boxes containing an annotated glomerulus in the center. Negative examples 
are 200 x 200 sub-images picked from random locations in the kidney microscopy images. Erom each sample, a 
512-dimensional R-HOG descriptor is extracted. 
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Training Data for Segmentation Stage. As described in the Segmentation subsection, the boundary likeliness is 
computed in every position on the m line segments. This boundary likeliness is the SVM score. The position lying 
on the true boundary of a glomerulus is considered as a positive example for the SVM, and the other positions are 
negative examples. To construct the training data for segmentation, the positive sub-images in the training dataset 
for pre-screening are reused. 

Training Data for Classification Stage. Examples in the training data for pre-screening are used again for 
training in the classification stage, but with a different set of features extracted via S-HOG. For each training data 
sample, the previously described segmentation algorithm estimates the boundary of the glomerulus. Based on the 
estimated boundary, the statistics of oriented gradients are computed to obtain S-HOG feature vectors. This procedure 
is done for both positive and negative examples, even though negative examples do not contain a glomerulus. 

Materials 

The images used in the present study had been generated in a previous study [S], and only an overview is given in 
this subsection. 

Male 6-week-old SD and SDT rats were purchased from GLEA, Inc. (Tokyo, Japan) and were housed with a 12-h 
light-dark cycle and free access to water and chow. 

Five SD and SDT Rats at 16 and 24 weeks of age were euthanized under ether anesthesia. Their kidneys were 
removed and immediately fixed in 10% neutralized buffered formalin. The formalin-fixed kidneys were embedded in 
paraffin. For immunohistochemistry, kidney paraffin sections were deparafhnized and incubated overnight at 4°G with 
anti-desmin mouse monoclonal antibody (Dako, Glostrup, Denmark) followed by horseradish peroxidase-conjugated 
secondary antibody (anti-mouse immunoglobulin goat polyclonal antibody; Nichirei, Tokyo, Japan). The sections were 
colorized brown with 3,3’-diaminobenzidine. Whole slide images of the sections were obtained with Aperio Scan Scope 
XT (Leica Microsystems, Wetzlar, Germany). All animal experiments were performed in accordance with the Act on 
Welfare and Management of Animals and the institutional guidelines, and approved by the institutional Gommittee 
of Animal Experiments of New Drug Development Research Genter Inc. (Hokkaido, Japan). 

3 Results and Discussion 

In this section, the detection performance is demonstrated by showing the experimental comparisons between S-HOG 
and R-HOG [Sl [9]. 

A set of 20 whole-kidney-section images is used in the experiments. The dataset is the same as the one used by 
[5]. The image sizes are 9,849 x 10,944 pixels in average. Each image is from one of four groups: 16-week-age SD 
rat, 16-week-age SDT rat, 24-week-age SD rat, and 24-week-age SDT rat. Henceforth, for simplicity, we will refer to 
them as 16SD, 16SDT, 24SD, 24SDT, respectively, each group containing five images. For performance evaluation, 
we annotated every glomerulus in the images manually. We divided the image set into five subsets: Set A, Set B, Set 
G, Set D, and Set E. Each subset consists of a 16SD image, a 16SDT image, a 24SD image, and a 24SDT image. For 
assessment of detection performance, the position of every glomerulus in the images is annotated and, for evaluation of 
segmentation performance, the areas of glomeruli in Set A and Set B are located manually using a graphics software. 

As described in the previous section, our method has three stages: pre-screening, segmentation, and classification. 
Each stage uses its own SVM trained with a hyper-parameter C. In the classification stage, a threshold 0 is used 
to classify an example; if the SVM score is over the threshold 6>, the example is predicted as positive, otherwise, 
negative. For the pre-screening and classification stages. Set A was used for training SVM, and Set B was exploited 
for determining the optimal combination of {C^O). Sets G, D, and E were for performance evaluation. SVM for the 
segmentation stage provides us with the boundary likelihood function. The regularization parameter C for the SVM 
is determined via the holdout method within Set A. Seventy percent of the glomeruli in Set A are randomly selected 
for training, and the rest is used for validation. The resulting parameter values were (C, 0) = (10, 2) for pre-screening, 
C = 10 for segmentation, and {C^O) = (10, —1.5) for classification. 

Detection Performance 

Figure 0] illustrates examples of detected glomeruli. In the two images, the candidate glomeruli passed through pre¬ 
screening are depicted with rings that represent the boundaries estimated in the segmentation stage. The numbers 
printed above the rings are the scores produced by SVM in the classification stage. Gandidate glomeruli with SVM 
scores below 0 = —1.5 are excluded from the final detection results. The excluded candidates are depicted with 
blue rings, and the remaining glomeruli with red rings. It can be observed that non-glomerulus areas are excluded 
effectively, whereas true glomeruli are estimated correctly. 


6 


For quantitative assessment of detection performance, true positives, false positives, and false negatives have to be 
defined. True positive glomeruli (TPG) are identified as correctly detected glomeruli, false positive glomeruli (FPG) 
are wrongly detected glomeruli, and false negative glomeruli (FNG) are the ones that could not be detected. From 
the definitions of TPG, FPG, and FNG, we can compute for the three widely used performance measures: F-measure, 
Precision, and Recall. Precision is the ratio of TPG to detected glomeruli (i.e. TPG/(TPG + FPG)), Recall is the 
ratio of TPG to true glomeruli (i.e. TPG/(TPG + FNG)), and F-measure is the harmonic mean of the Precision and 
the Recall. 

Figure [S] shows the plots of the F-measure, Precision, and Recall for each testing image. S-HOG achieves an average 
of 0.866, 0.874, and 0.897 for F-measure, Precision, and Recall, respectively, whereas R-HOG obtained 0.838, 0.777, 
and 0.911, respectively. In applying detection methods to pathological evaluation. Precision is more important than 
Recall [H], and in this study, S-HOG achieved considerably higher Precision at a small sacrifice of Recall. Two-sample 
t-test is performed to assess the differences statistically. While no statistical difference of Recall can be detected 
(P-value is 3.47 • 10“^), the differences among F-measure and Precision are significant (P-values is 1.34 • 10“^ and 
3.75-10“^, respectively). 

Segmentation Performance 

One advantage of the proposed method is that the boundaries of the detected glomeruli can be obtained. These 
boundaries provide useful information for pathological evaluation [9]. Here we give a discussion of the performance 
of the segmentation algorithm. To quantify the accuracy of the estimated areas within the predicted boundaries, 
993 annotated glomeruli in Set B were used. True positive area (TPA), false positive area (FPA), and false negative 
area (FNA) are defined as follows: TPA is the intersection of the true area and estimated area; FPA is the relative 
complement of the true area in the estimated area; FNA is the relative complement of the estimated area in the true 
area. For each glomerulus and its estimated area, F-measure, Precision, and Recall can be obtained by counting the 
pixels in the TPA, FPA, and FNA. The histograms of the F-measure, Precision, and Recall are plotted in Figure [6l 
where the frequency is normalized so that the integral is one. Among the glomeruli, 90.1% are estimated to have 
F-measure over 0.8, ensuring reliable assessment of medicinal effect for pharmaceutical development. 

The computational time of the new segmentation algorithm, DCDP, is compared with that of EDP. Note that the 
two algorithms solve the same optimization problem, and it can be shown that both algorithms always find an exact 
optimal solution. DCDP and EDP are implemented in C++ language, and the runtimes are measured on a Linux 
machine with Intel(R) Core(TM) i7 CPU and 8Gb memory. Eirst, the number of times when the O(nm^) DP routine 
was invoked, which we denote by ridp, is counted using the annotated glomeruli in Set B. Eigure[7^ shows the box-plot 
of Udp for all methods. While the value of Udp for EDP is always n, the values for DCDP depend on the input images 
and the splitting schemes. Half Split, Max Split, and Adap Split (Section [B]). Eor 46.32% of glomeruli, the optimal 
solutions are found within the first DP routine (i.e. Udp = !)• The medians of the n^pS when using Half Split, Max 
Split, and Adap Split are 5, 3, 3, respectively. In other words, the medians of the depths of the branching tree for each 
scheme, respectively, are 3, 2 and 2. The 75 percentiles of Udp’s are 11, 7, 5, respectively. Eor no glomeruli glomeruli, 
ridp of Adap Split is larger than n, whereas the number of glomeruli with n^p > n are 4 (0.40%) and 16 (1.61%) for 
Half Split and Max Split. This implies that Adap Split is the smartest heuristic among the three splitting schemes. 
As considered in Subsection m Adap Split produces the same solution pl in the branches less frequently than the 
three other schemes. In Half Split and Max Split, the frequencies (# of glomeluli) for the cases that the solution pl 
in the top branch appears again in the second branches is 414 and 314, respectively. Those numbers are much larger 
than the frequency in Adap Split which is only 97. This explains why Adap Split is faster. The actual runtimes of 
each methods are depicted in Eigure [7)3, where the medians of the computational times are 0.0866, 0.0570, 0.0560, 
and 0.418 msec, respectively, for Half Split, Max Split, Adap Split and EDP. The 75 percentiles of the computational 
times are 0.171, 0.117, 0.0856, and 0.426 msec, respectively. Since these values are proportional to the Udp’s, then the 
ratios among the runtimes are almost same as the ratios among the n^pS. These results conclude that the proposed 
algorithm DCDP achieves an exact optimal solution much more efficiently than the existing algorithm EDP solves the 
same problem, and Adap Scheme is the fastest splitting scheme. 

4 Conclusions 

In this paper, a new descriptor. Segmental HOG, was proposed for specific organ detection in microscopy images. The 
descriptor was based on the boundary of glomeruli to acquire robustness to variations of intensities, sizes, and shapes. 
A new segmentation algorithm, DCDP, was developed to locate the boundary of possible candidates of glomeruli. 
Empirical results show significant improvement compared to the state-of-the-art descriptor. Rectangular HOG, for the 
task of glomerulus detection in microscopy images. Moreover, experimental results reveal that DCDP is much faster 
than the existing segmentation algorithm EDP. 
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Several possible extensions of the proposed method can be considered. For instance, appropriate size of the sliding 
window should be chosen if the proposed method is applied to microscopic images with different resolutions. Also, 
while the boundary likeliness function is the same for any direction in the segmentation algorithm, different boundary 
likeliness functions can be used for detecting other organs that have orientation. As for the block division of the S-HOG 
descriptor, 24 blocks are used in this study as depicted in Figure [3t, but a different number of blocks with a different 
division can be used for another application. Future work includes exploring such extensions in other applications of 
Segmental HOG. 


A Pruning 


Pruning can accelerate the DGDP algorithm. Gonsider the case where the lower bound such that 


max 

pes{Nn) 


j{p) > e, 


is known in advance when searching for the solution in S{Xo). If i > maXp^^^^^Xo) J{p)^ then no optimal solution is in 
S{Xq) since 


max J(p) > £ > 

pesiNn) 


max J(p) > max J(p). 

peSi^iXo) pesiXo) 


Based on this fact, the pruning step is added to obtain Algorithmic and we have the following lemma. 


Lemma 1. For any subset Xq C Nn and G MU {—oo}, when the algorithm runs with (po, = DCDP(Xo,^), 

the returned tuple (po, satisfies one of the following: 

• Case G: If maXp^^^Xo) J{p) > £^ then 


Po G argmax J(p), and Jq = £o = J(po) > £• 
peSiXo) 


• Case L: If maXp^^^Xo) J{p) < then Jq < Iq = £. 


Algorithm A (po, Jo,£o) = DCDP(Xo,^). Modification of DCDP.Basic by adding pruning steps. 
Input: Xq C Nn- 
1: PL e argmaxpe5^(2:^) J(p); 

2: if ^ > J(pl) then 
3: {Rule A} 

4: Jo := —oo; £o := £; 

5: return; 

6: end if 

7: if pl G S{Xo) then 
8: {Rule B} 

9: Po :=Pl; Jo ■= J{pl)-, 4 := max(^, Jo)); 

10: return; 

11: end if 
12: {Rule C} 

13: (Xi, J 2 ) := Split(Xo); 

14: (PI, Ji,4) :=DCDP(Xi,0; 

15: (P2,J2,^2) := DCDP(X2,i?i); 

16: i* e argmax^gii 2} Ji'^ 

17: Po ;=Pi*; Jo ;= J(Pi*); 

18 : £0 '■= max(^ 2 , Jo)', 


Section [O gives the proof of the lemma. Finally, from Lemma [H we can derive the following theorem which is an 
important theoretical result of this study. 

Theorem 1. The optimal solution of the problem 0 is obtained by invoking (p*, = DCDP(Nn, —oo). 

B Splitting Schemes 

For the implementation of Split(Xo), we considered three schemes: Half Split, Max Split, and Adap Split. 
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Half Split. In this scheme, the subset of indices Xq is simply divided into the first half and the second half. The 
resulting Xi and X 2 are sets of consecutive integers. For instance, this scheme divides Xq = {7,8,9,10,11,12} into 
Xi = {7, 8, 9} and X 2 = {10,11,12}. To increase the lower-bound £ defined earlier, a heuristic that swaps Xi with X 2 if 

maxX^(h) < maxX^(h) 


is employed. 

Max Split. Similar to Half Split, Xi and X 2 generated by Max Split are sets of consecutive integers, although the 
splitting points are different. The splitting point of Half Split is the center of the interval Xq, whereas the splitting 
point of Max Split is given by h*(Xo) := argmax^^j^ Lm{h). For example, if h^(Xo) = 8 for Xq = {7,8,9,10,11,12}, 
this scheme outputs Xi = {7, 8} and X 2 = {9,10,11,12}. In general, the resulting divisions are given by 

Xi := {h G Xo I h < K{Xq)} , and X 2 := {h G Xq | h > K{Xq)} . 

If card(X 2 ) = 0, the entry h^{Xo) is moved from Xi to X 2 . A heuristic that swaps Xi with X 2 if card(Xi) > card(X 2 ) is 
applied. 

Adap Split. Different from the above two schemes, Adap Scheme adaptively determines the splitting point using 
the current solution pl,o •= ^^S^^^pg5l(Xo) denote the first entry and the last entry in the vector pl,o 

by Pi and respectively. The splitting point of Adap Scheme is given by the center between Pi and p^. For 
instance, when p? = 9 and p^ = 12 Xq = {7, 8, 9,10,11,12}, this splitting scheme divides Xq into Xi = {7, 8, 9,10} and 
X 2 = {11,12}. In general, the resulting divisions are given by 

Xi := |/1 e lo I< — 

The smallest entry in X 2 is moved to Xi if card(Xi) = 0. The largest entry in Xi is moved to X 2 if card(X 2 ) =0. A 
swapping heuristic used in Max Split is then applied. 

B.l Why Adap Split is better 

Adap Scheme is expected to be the smartest heuristic among the three splitting schemes. To describe the reason, let 
us illustrate the process of DCDP on a small toy problem with (n, m,<^) = (12,8,1) shown in Figure [HI The original 
problem and the relaxed problem are depicted in Figure |8}i,b, respectively. When running DCDP(N8, — 00 ), it is 
observed that pl,o •= ^(p) ^ ‘^(Xo), and thereby the set Xq is divided into Xi and X 2 to produce two 

new branches DCDP(Xi, — oc) and DCDP(X 2 ,-^i) where ii will be computed by the former branch DCDP(Xi, — 00 ). 

If using the Adap Split scheme, the two subsets are Xi = {7,8} and X 2 = {1,2,3,4, 5,6}. In the branch of 
DCDP(Xi, — 00 ), pl,i •= J{p) shown in Figure[8t), implying that pl,i is the maximizer 

of J{p) over 5(Xi) and no more branching occurs. ii = J(pl,i) = 10.1 is set. Next, DCDP(X 2 ,-^i) is invoked and 
Pl ,2 •= ^(P) computed in the branch (Figure [H]l). It is then observed that J(pl, 2 ) = 9.5 < 10.1 = -£ 1 , 

which implies that the optimal solution is not in S{X 2 ). Hence, pl,i turns out the optimal solution of the original 
solution. 

Meanwhile, if Half Split is applied, the set Xq = Ns is divided into Xi = {5, 6, 7, 8} and X 2 = {1, 2, 3,4} (Figure[8t,f). 
In the branch of DCDP(Xi, — 00 ), the solution of the relaxed problem is again pl,i •= J{p) — Pl,o 0 

5(Xi), leading to further branching of this branch (Figure |8t). Actually, in our experiments described in Results 
section, it is observed that Half Split and Max Split frequently encounter pl,i = Pl,o 01 * Pl,2 = Pl,o- Whenever 
Pl,i = Pl,0 5 further new branches for the divisions of Xi are produced, because pl,i = Pl,o ^ S{Xq) leading to 
Pl,i ^ 5 (Xi) C tS(Xo). Similarly, when pl,2 = Pxpi further new branches for the divisions of Xi are always born. 

In Adap Split, as desribed in Results section, pl,i = Pl,o ot pl ,2 = Pl,o is less likely to happen, resulting in less 
branches. 


q. 


and I 2 := <h £lo\h > 


P°i 




C Proof of Lemma [T] 


We shall use the following notation: For any X C Nn, 


J(X) := max J{p), and Jl{^) •= ffiax J{p). 
pes(x) pGvSl(x) 
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max(J(Ii), J(l 2 )), 


The following relationships will be used in this proof: 

J(Pl) ^ Jl(Xo) J(Xo) ^ 

where the labels (eqA), (ineqA), and (eqB) are used to distinguish these equalities and the inequality in later descrip¬ 
tions of this proof. The inequality follows from the fact that S(Xq) C 5l(To), while the second equality eqB follows 
from5(Xo) =5 (Xi)U5(X2). 

We will prove the lemma by induction. For the case where card(Xo) = 1, observe that 5 l(To) = S(Xq)^ implying 
that 


Pl e argmax J(p) and J(Xo) = Jl(53) = Jq- 
pes{Xo) 

If J(Xo) = J(pl) < '^5 then by Rule A in the algorithm, Jq = — oo and Jo < On the other hand, if J(pl) ^ 

then since pl ^ have pl ^ ‘5(Xo)- Thus, by Rule B, we have po = Pl and Jo = ^ 

Therefore, the lemma is true for card(Xo) = 1- 

Let us now assume that the lemma holds for any Xq C Nn such that card(Xo) < k to show that the lemma is also 
established for any Xq such that card(Xo) = k. Now suppose that Xq C Nn and card(Xo) = k. The following is an 
exhaustive list of all possible cases: 

(1) £ > Jl(Xo), 

(2) Jl(Jo) > £ and pl G 5(Xo), 

(3) Jl(Xo) > J(Xo) = J(Xi) = J(X 2 ) > £ and pl ^ 5(Xo), 

(4) Jl(Xo) > J(Xo) = J{Ti) > £ > J{X 2 ) and pl ^ 5(Xo), 

(5) Jl(To) > J(Xo) = J(Xi) > J{X 2 ) > £ and pl ^ 5(Xo), 

(6) Jl(To) > J(Xo) = J{T 2 ) > J{Xi) > £ and pl ^ 5(Xo), 

(7) MTo) > J{Xo) = J{X 2 ) > £ > JiXi) and pl ^ 5(Xo), 

(8) Jl(Jo) > £ > J{Xo) and pl 0 S{Xo). 

We shall show that for each of the seven cases above, either Case G or Case L is true. 

/ A") (ineqA) 

(1) Since £ > Jl(Jo) ^ ^ J(pl) > then by Rule A, we have Jo = —oo and Jq < £o = £• Thus, Case L is 

satisfied. 

(2) If Jl(Jo) > £, then J(Pl) > £. Moreover, pl G S{Xo) implies pl G argmax^^^j^j^) J(p) and J(pl) = J(Jo). 

Then by Rule B, Jo = J(pl) = J(Jo) = £o and Po = Pl G argmaXp^^^^j^^ J{p)‘ Therefore, Case G holds. 

i 0Q 

(3) Given J(pl) = Jl(Jo) > £ and pl ^ S{Xq)^ Rule C is applied. Observe that (pi, Ji,^i) satisfies Case G since 

J(Xi) = maXp^^^Xi) J{p)^ £‘ Thus, Ji = £i = J{Xi) > £. Similarly, (p 2 , J 2 ,'^ 2 ) also satisfies Case G, and we 
have J 2 = £2 = J{^ 2 ) ^ J{^i) = '7(Xo)- Therefore, the returned value po of DCDP(Xo 5 '^) is equal to pi or to 

P 2 , both of which are in argmaXp^^^'j^^ J{p)' Furthermore, we obtain £0 = Jq = Ji = J 2 = J(Jo) > £• Hence, 
we have Case G. 

(4) As with Case (3), (pi, Ji^£i) satisfies Case G, and we have pi G argmaXp^^^^j^^ J(p) and Ji = £i = J{Xi) > £. 
An so it follows that maXp^^^Xs) J{P 2 ) = J{^ 2 ) < £ < £\^ and Case L holds for (p 2 , J 2 A 2 ) with J 2 < £2 = £\. 
Therefore, output po = Pi, which is in argmaXp^^^^Xo) '^(p) since iS(Xi) C ^(Xo), and Jo = Ji = £2 = £0 ^ £• 
This gives us Case G. 

(5) In a similar logic as in Case (4), Case G is true for (pi, Ji,^i), and pi G argmaXp^^^^Xi) J{p) and Ji = £1 = 
J{^i) ^ Given that J{X 2 ) < J(Xi), then J(X 2 ) < ^ 1 , hence. Case L also holds for (p 2 , J 2 ,'^ 2 ) in this sub-case. 
Therefore, the same conclusion from Case (4) follows. 

(6) Since px ^ S{Xo)^ Rule C is implemented. Note that with maxp^^^x^) J(p) = J(Xi) > (pi, Ji,^i) follows 
Case G, and we have Ji = £i = J(Xi) > £. This implies that J(X 2 ) > J(Xi) = £ 1 . Therefore, Case G also 
holds for (p2, J2,'^2), and we obtain J 2 = £2 = J(X2) > J(Xi) = Ji. Thus, DCDP(Xo 5 '^) outputs Po = P2 G 
argmaXp^^^'Xo) ^ip) ^^4 Jq = J 2 = £2 = £0 > £• Hence, Case G holds. 
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(7) Similar to the previous cases, we apply Rule C. Given £ > J{Ti) = maXp^^^Xi) J{p)i by Case L, Ji < £i = i. 

and so it follows that J{T 2 ) > ii = i. Therefore, (p 2 , -^ 2 ,-^ 2 ) satisfies Case G, with J 2 = £2 = J{^ 2 ) > £• Hence, 
the returned values of DCDP(Xo,-^) are Po = P2 ^ ^(p) Jo = J 2 = £2 = £0 > £^ and Case G 

is satisfied. 

(8) Likewise, we implement Rule C for this case. Observe that since max( J(Xi), J(X 2 )) J{Jo) < £, then both 
(pi, Ji^£i) and (p 2 , J 2 ^£ 2 ) satisfy Case L. Thus, we have Ji < £i = £ and J 2 < £1 = £ 2 - Therefore, DCDP(Xo,^) 
outputs Jo = max(Ji, J 2 ) < £2 = £0 = h = £^ which corresponds to Case L. 
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(a) Pre-screening by slid- (b) Extraction of Seg- 
ing window. mental HOG features and 

classification 


Figure 1: Flow of Our Method. In this study, a new descriptor, Segmental HOG, is developed to detect glomeruli 
in huge microscopy images. SVM is combined with Segmental HOG to classify candidates of glomeruli that passed 
the pre-screening stage. 


(a) 



(b) 



Figure 2: Candidate Glomerulus and Line Segments. S-HOG is based on the boundary of the objects of interest. 
Suppose the boundary of a candidate glomerulus (Panel(a)) is to be located. Boundary likeliness is computed at every 
point on m{= 36) line segments placed uniformly in all m directions (Panel(b)). The boundary likeliness is computed 
at n points on each line segment. The n locations are depicted with dots in Panel (b). 
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(a) Blocks of R-HOG 



(b) Feature Vector of R-HOG 








(c) Blocks of S-HOG 



(d) Feature Vector of S-HOG 


Figure 3: Rectangular HOG and Segmental HOG. R-HOG has been used for object detection in a wide variety 
of applications. R-HOG is the concatenation of statistics in each block in a grid dividing a rectangular region. On the 
other hand, the blocks of the proposed descriptor, S-HOG, are based on the segmentation of the object of interest. 




Figure 4: Examples of Detected Glomeruli. The estimated boundaries of the glomeruli are depicted with red rings. 
The areas surrounded with blue rings are passed through the pre-screening stage, but removed in the classification 
stage. The numbers are the SVM scores resulting from the classification stage. The areas with SVM scores over 
0 = —1.5 are classified as a glomerulus. It can be observed that false positives such as vessels detected in the 
pre-screening stage were successfully removed in the classification stage. 
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Figure 5: Detection Performances. The proposed descriptor, S-HOG, achieves evident improvement in F-measure 
compared to the existing descriptor, R-HOG. With small loss of true positives, S-HOG halves false positives of R-HOG. 
(See subsection on Detection Performance for details.) 



F-measure Precision Recall 


Figure 6: Segmentation Performances. The number of glomeruli are tallied to make a histogram with the F- 
measure. Precision, and Recall of the pixels based on comparison of true segmentation with estimated segmentation 
on the x-axes. (See Segmentation Performance subsection for details.) 


(a) (b) 




Figure 7: Runtime Comparisons. In Panel(a), ndp of DGDP with three splitting schemes and EDP is shown, 
where ndp is the number of invoking the O(mn^) DP routine. The computational time of each algorithm is plotted in 
Panel(b). 
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(a) 5(N8) 



(b) 5 l(N8) 



(c).Sl({7,8}) 



(d)5L({l,2,...,6}) 



(e) 5l({5,6,7,8}) 



(f)5L({l,2,3,4}) 



Figure 8: Relaxed Problems in DCDP. Here a segmentation problem dU with (n, m,^) = (12,8,1) is considered. 
In Panel (a), the sizes of red circles indicate the quantities of boundary likeliness. The feasible region S{Xq) is the 
set of polygons, where Xq = N^. Overlapping all the polygons yields the gray edges. The optimal polygon is drawn 
with the red edges. DCDP relaxes the feasible region to get Si,{Xo). In Panel (b), the relaxed feasible region 5 l(Xo) is 
depicted. The blue polygon is the optimal solution for the relaxed problem pl,o = ^(p)* relaxed 

problems of four sub-problems with 5l({7, 8}), 5l({ 1, 2,..., 6}), iSL({5, 6, 7, 8}), and 5l({1, 2, 3,4}) are illustrated 
in Panels (c),(d),(e),(f). The blue polygons in (c),(d),(e),(f) are the optimal solution of the four relaxed problems, 
respectively. 
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